##################################################
# Replication Code
# Taeyong Park and Andrew Reeves
# "Local Unemployment and Voting for President: Uncovering Causal Mechanisms"
# Summary: Code for running sensitivity analysis and replicating Figure 1 in the online appendix.
##################################################

rm(list = ls())
library(foreign); library(mediation)

##########
## 2008 ##
##########
data08 = read.dta("data08Imputed1.dta") 

# Evaluations
data08$VBad = rep(0, nrow(data08))
data08$VBad[data08$natecon5==5] = 1
data08$Bad = rep(0, nrow(data08))
data08$Bad[data08$natecon5==4] = 1
data08$Good = rep(0, nrow(data08))
data08$Good[data08$natecon5==2] = 1
data08$VGood = rep(0, nrow(data08))
data08$VGood[data08$natecon5==1] = 1

# Age
data08$younger30 = rep(0, nrow(data08))
data08$younger30[data08$age<30] = 1
data08$older65 = rep(0, nrow(data08))
data08$older65[data08$age>=65] = 1

# Race
data08$black = rep(0, nrow(data08))
data08$black[data08$race==1]=1
data08$latino = rep(0, nrow(data08))
data08$latino[data08$race==2]=1

# Employment status
data08$fullTime = rep(0, nrow(data08))
data08$fullTime[data08$employment == 1] = 1 
data08$partTime = rep(0, nrow(data08))
data08$partTime[data08$employment == 2] = 1 
data08$unemployed = rep(0, nrow(data08))
data08$unemployed[data08$employment == 3] = 1 

# Education
data08$noHighSchool = rep(0, nrow(data08))
data08$noHighSchool[data08$educ==4] = 1
data08$HighSchool = rep(0, nrow(data08))
data08$HighSchool[data08$educ==3] = 1
data08$someCollege = rep(0, nrow(data08))
data08$someCollege[data08$educ==1] = 1
data08$fourCollege = rep(0, nrow(data08))
data08$fourCollege[data08$educ==2] = 1

# Party ID
data08$Rep = rep(0, nrow(data08))
data08$Rep[data08$party3==3] = 1
data08$Ind = rep(0, nrow(data08))
data08$Ind[data08$party3==2] = 1
data08$Dem = rep(0, nrow(data08))
data08$Dem[data08$party3==1] = 1

# News interest
data08$Very=ifelse(data08$newsInt==4, 1, 0)
data08$Some=ifelse(data08$newsInt==3, 1 ,0)
data08$Few=ifelse(data08$newsInt==1 | data08$newsInt==2, 1, 0)

# Unemployment
data08$avgUnemp07 = apply(data08[,7:18], 
                        FUN=mean, MARGIN=1)
data08$avgUnemp08 = apply(data08[,19:30], 
                        FUN=mean, MARGIN=1)
data08$chgUnemp0807=data08$avgUnemp08-data08$avgUnemp07

# Income 
data08$chgInc0807=(data08$inc08-data08$inc07)/1000
data08$inc08=data08$inc08/1000

# Gas price
data08$avgGas = (data08$gasOct08 + data08$gasSept08 + data08$gasAug08)/3

# Foreclosure
data08$totForcAugOct = data08$forcAug08  + data08$forcSep08+   data08$forcOct08

modelM08 = lm(natecon5 ~ chgUnemp0807 +  avgUnemp08 + chgInc0807 + inc08 + 
                 avgGas + totForcAugOct +
                 younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                 income + unemployed + fullTime + partTime + ownHome +
                 Very + Few + ideol + Dem + Rep + Ind, data=data08)
modelY08 = glm(pvote2 ~ natecon5 + chgUnemp0807 +  avgUnemp08 + chgInc0807 + inc08 + 
                  avgGas + totForcAugOct +
                  younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                  income + unemployed + fullTime + partTime + ownHome +
                  Very + Few + ideol + Dem + Rep + Ind, data=data08, family=binomial("probit"))

mediation08 = mediate(modelM08, modelY08, treat="chgUnemp0807", mediator="natecon5", sims=500)
sens08 = medsens(mediation08, rho.by=.1)



##########
## 2012 ##
##########
data12 = read.dta("data12Imputed1.dta") 

# Evaluations
data12$VBad = rep(0, nrow(data12))
data12$VBad[data12$natecon5==5] = 1
data12$Bad = rep(0, nrow(data12))
data12$Bad[data12$natecon5==4] = 1
data12$Good = rep(0, nrow(data12))
data12$Good[data12$natecon5==2] = 1
data12$VGood = rep(0, nrow(data12))
data12$VGood[data12$natecon5==1] = 1

# Age
data12$younger30 = rep(0, nrow(data12))
data12$younger30[data12$age<30] = 1
data12$older65 = rep(0, nrow(data12))
data12$older65[data12$age>=65] = 1

# Race
data12$black = rep(0, nrow(data12))
data12$black[data12$race==1]=1
data12$latino = rep(0, nrow(data12))
data12$latino[data12$race==2]=1

# Employment status
data12$fullTime = rep(0, nrow(data12))
data12$fullTime[data12$employment == 1] = 1 
data12$partTime = rep(0, nrow(data12))
data12$partTime[data12$employment == 2] = 1 
data12$unemployed = rep(0, nrow(data12))
data12$unemployed[data12$employment == 3] = 1 

# Education
data12$noHighSchool = rep(0, nrow(data12))
data12$noHighSchool[data12$educ==4] = 1
data12$HighSchool = rep(0, nrow(data12))
data12$HighSchool[data12$educ==3] = 1
data12$someCollege = rep(0, nrow(data12))
data12$someCollege[data12$educ==1] = 1
data12$fourCollege = rep(0, nrow(data12))
data12$fourCollege[data12$educ==2] = 1

# Party ID
data12$Rep = rep(0, nrow(data12))
data12$Rep[data12$party3==3] = 1
data12$Ind = rep(0, nrow(data12))
data12$Ind[data12$party3==2] = 1
data12$Dem = rep(0, nrow(data12))
data12$Dem[data12$party3==1] = 1

# News interest
data12$Very=ifelse(data12$newsInt==4, 1, 0)
data12$Some=ifelse(data12$newsInt==3, 1 ,0)
data12$Few=ifelse(data12$newsInt==1 | data12$newsInt==2, 1, 0)


# Unemployment
data12$avgUnemp11 = apply(data12[,15:26], 
                        FUN=mean, MARGIN=1)
data12$avgUnemp12 = apply(data12[,27:38], 
                        FUN=mean, MARGIN=1)
data12$chgUnemp1211=data12$avgUnemp12-data12$avgUnemp11

# Income
data12$chgInc1211=(data12$inc12-data12$inc11)/1000
data12$inc12 = data12$inc12/1000

# Gas price
data12$avgGas = (data12$gasOct + data12$gasSep + data12$gasAug)/3

modelM12 = lm(natecon5 ~ chgUnemp1211 +  avgUnemp12 + chgInc1211 + inc12 +
                 avgGas + totForcAugOct +
                 younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                 income + unemployed + fullTime + partTime + ownHome +
                 Very + Few + ideol + Dem + Rep + Ind, data=data12)
modelY12 = glm(pvote2 ~ natecon5 + chgUnemp1211 +  avgUnemp12 + chgInc1211 + inc12 + 
                  avgGas + totForcAugOct +
                  younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                  income + unemployed + fullTime + partTime + ownHome +
                  Very + Few + ideol + Dem + Rep + Ind, data=data12, family=binomial("probit"))

mediation12 = mediate(modelM12, modelY12, treat="chgUnemp1211", mediator="natecon5", sims=500)
sens12 = medsens(mediation12, rho.by=.1)



##########
## 2016 ##
##########
data16 = read.dta("data16Imputed1.dta") # Import one imputed dataset. Must check if the analysis results are robust to other datasets

# Evaluations
table(data16$natecon5)
data16$VBad = rep(0, nrow(data16))
data16$VBad[data16$natecon5==5] = 1
table(data16$VBad)

data16$Bad = rep(0, nrow(data16))
data16$Bad[data16$natecon5==4] = 1
table(data16$Bad)

data16$Good = rep(0, nrow(data16))
data16$Good[data16$natecon5==2] = 1
table(data16$Good)

data16$VGood = rep(0, nrow(data16))
data16$VGood[data16$natecon5==1] = 1
table(data16$VGood)

# Age
table(data16$age)
data16$younger30 = rep(0, nrow(data16))
data16$younger30[data16$age<30] = 1
table(data16$younger30)

data16$older65 = rep(0, nrow(data16))
data16$older65[data16$age>=65] = 1
table(data16$older65)

# Race
table(data16$raceNew)
data16$black = rep(0, nrow(data16))
data16$black[data16$raceNew==1]=1
table(data16$black)
data16$latino = rep(0, nrow(data16))
data16$latino[data16$raceNew==2]=1
table(data16$latino)

# Employment status
table(data16$employment)
data16$fullTime = rep(0, nrow(data16))
data16$fullTime[data16$employment == 1] = 1 
table(data16$fullTime)
data16$partTime = rep(0, nrow(data16))
data16$partTime[data16$employment == 2] = 1 
table(data16$partTime)
data16$unemployed = rep(0, nrow(data16))
data16$unemployed[data16$employment == 3] = 1 
table(data16$unemployed)

# Education
table(data16$educNew)
data16$noHighSchool = rep(0, nrow(data16))
data16$noHighSchool[data16$educNew==4] = 1
table(data16$noHighSchool)

data16$HighSchool = rep(0, nrow(data16))
data16$HighSchool[data16$educNew==3] = 1
table(data16$HighSchool)

data16$someCollege = rep(0, nrow(data16))
data16$someCollege[data16$educNew==1] = 1
table(data16$someCollege)

data16$fourCollege = rep(0, nrow(data16))
data16$fourCollege[data16$educNew==2] = 1
table(data16$fourCollege)

# Party ID
table(data16$party3)
data16$Rep = rep(0, nrow(data16))
data16$Rep[data16$party3==3] = 1
table(data16$Rep) 
data16$Ind = rep(0, nrow(data16))
data16$Ind[data16$party3==2] = 1
table(data16$Ind) 
data16$Dem = rep(0, nrow(data16))
data16$Dem[data16$party3==1] = 1
table(data16$Dem) 


# News interest
table(data16$newsInt)
data16$Very=ifelse(data16$newsInt==4, 1, 0)
data16$Some=ifelse(data16$newsInt==3, 1 ,0)
data16$Few=ifelse(data16$newsInt==1 | data16$newsInt==2, 1, 0)
table(data16$Very)
table(data16$Some)
table(data16$Few)

# Unemp
data16$avgUnemp15 = apply(data16[,5:16], 
                        FUN=mean, MARGIN=1)
data16$avgUnemp16 = apply(data16[,17:28], 
                        FUN=mean, MARGIN=1)
data16$chgUnemp1615=data16$avgUnemp16-data16$avgUnemp15

# Income
data16$chgInc1615=(data16$inc16-data16$inc15)/1000
data16$inc16 = data16$inc16/1000


modelM16 = lm(natecon5 ~ chgUnemp1615 + avgUnemp16 + chgInc1615 +  inc16 + 
                 younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +
                 income + ownHome + newsInt + ideol + Dem + Rep + Ind, data=data16)
modelY16 = glm(pvote2 ~ natecon5 + chgUnemp1615 + avgUnemp16 + chgInc1615 +  inc16 +
                  younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +
                  income + ownHome + newsInt + ideol + Dem + Rep + Ind, data=data16, family=binomial("probit"))

mediation16 = mediate(modelM16, modelY16, treat="chgUnemp1615", mediator="natecon5", sims=500)
sens16 = medsens(mediation16, rho.by=.1)


##########
## PLOT ##
##########
plot(sens08, main="2008", ylim=c(-0.01, 0.01))
plot(sens08, sens.par="R2", r.type="residual", sign.prod="negative", main="2008")

plot(sens12, main="2012", ylim=c(-0.03, 0.03))
plot(sens12, sens.par="R2", r.type="residual", sign.prod="negative", main="2012")

plot(sens16, main="2016", ylim=c(-0.03, 0.03))
plot(sens16, sens.par="R2", r.type="residual", sign.prod="negative", main="2016")


